A fully relativistic lattice Boltzmann algorithm 
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Starting from the Maxwell- Juttner equilibrium distribution, we develop a relativistic lattice Boltz- 
mann (LB) algorithm capable of handling ultrarelativistic systems with flat, but expanding, space- 
times. The algorithm is validated through simulations of quark-gluon plasma, yielding excellent 
agreement with hydrodynamic simulations. The present scheme opens the possibility of transferring 
the recognized computational advantages of lattice kinetic theory to the context of both weakly and 
ultra-relativistic systems. 
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I. MOTIVATION 

The great success of the Relativistic Heavy-Ion Col- 
lider (RHIC) experimental program [lJ-Q has provided 
the motivation to come up with realistic and quantita- 
tive simulations of heavy-ion collisions. 

Since the bulk of particles produced in relativistic 
heavy- ion collisions is described by fluid dynamics Q, 
the center-piece of any complete simulation attempt will 
involve a viscous fluid dynamics algorithm. The major- 
ity of presently available fluid dynamics codes is able to 
handle smooth initial conditions in 2+1 dimensions in 
the presence of shear viscosity [Hfllj- However, it has by 
now been understood that the presence of event-by-event 
fluctuations in the initial state can lead to significantly 
different quantitative results with respect to smooth ini- 
tial conditions and may in some cases even explain 
qualitatively new phenomena. To be more specific, the 
presence of event-by-event fluctuations is the source of 
the non-vanishing elliptic flow found in RHIC experi- 
ments at central collisions, the source of hydrodynamic 
flow-fluctuations, and may (through the presence of so- 
called triangular flow V3) naturally explain the presence 
of the 'ridge phenomenon' found in experiments [L2Hl5j . 
Thus, it is probably fair to say that without including 
the effect of event-by-event fluctuations, a description 
of the medium created in heavy-ion collisions cannot be 
regarded as realistic. This provides the motivation to 
develop a fully relativistic and computationally efficient 
viscous fluid dynamics algorithm that can handle initial 
state fluctuations. Also, such an algorithm can be used 
to validate the only available 3+1 dimensional relativistic 
viscous hydrodynamics code by the McGill group [l6j . 
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Further motivation is provided by other systems where 
relativistic viscous fluid flows are of interest, such as as- 
trophysical systems and condensed matter systems such 
as graphene [17[. One particular question that arises in 
all this different systems is when relativistic fluid flow 
becomes turbulent, which involves a determination of 
the critical Reynolds number and the turbulent spectrum 
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ii. LATTICE KINETIC APPROACH TO 
HYDRODYNAMICS 

Fluid turbulence, both classical and relativistic, sets 
one of the most compelling challenges in modern com- 
putational physics. This motivates a relentless search 
for new and ever more efficient methods for solving the 
hydrodynamic equations of motion in the high-Reynolds 
turbulent regimes. In the last two decades, a new com- 
putational paradigm has emerged, which is based on the 
idea of solving hydrodynamic problems from the stand- 
point of Boltzmann kinetic theory. Apparently, this is 
rather counterintuitive, because the Boltzmann equation 
lives in a double-dimensional (phase) space, consisting 
of three dimensions in ordinary space, plus three addi- 
tional dimensions in momentum (velocity space). In ad- 
dition, the Boltzmann equation contains the details of 
microscopic interactions through a very complicated in- 
tegral collision operator in velocity space, which is com- 
putationally very demanding. As a result, the Boltz- 
mann equation has never been considered a practical 
tool for computational fluid dynamics, apart from the 
special case of rarefied gas dynamics, for which ordinary 
fluid dynamics is known to be inadequate. In the last 
two decades, however, minimal realizations of the Boltz- 
mann equation have been developed, which relinquish 
the aforementioned problems, and gave rise to a compu- 
tational method of remarkable elegance and outstanding 
computational efficiency. Since these minimal forms of 
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Boltzmann kinetic equations are formulated in a discrete 
velocity and space-time lattice, they have come to be 
known as Lattice Boltzmann Equation(s) (LBE) p8l - l30l | . 
To date, LB methods have met with amazing success 
across virtually all sectors of non-relativistic fluid dy- 
namics, from flows to porous media, to turbulent flows in 
complex geometries, multiphase, colloidal, hemodynamic 
flows, and magnetohydrodynamics 31-38]. However, rel- 
ativistic formulations have come into existence only very 
recently [3^, H(|. Indeed, the lattice formulation of the 
rclativistic Boltzmann equation presents a series of the- 
oretical and methodological challenges which have no 
counterpart in the non-relativistic realm. To date, some 
of these challenges have been bypassed by formulating 
a relativistic LBE (RLBE) top-down, i.e. by recover- 
ing the equations of relativistic hydrodynamics through 
a moment-matching procedure in the lattice. This gives 
rise to a very efficient computational scheme for mildly 
rclativistic systems, but does not guarantee the so-called 
"realizability" of RLBE, i.e. the fact that RLBE should 
be derived by an underlying microscopic model, or, at 
least, by a continuum version of a relativistic kinetic 
equation. Even leaving aside computational considera- 
tions, this is an important task in the process of placing 
the RLBE onto a solid conceptual framework. This is 
precisely the task accomplished in the present paper. 



III. RELATIVISTIC KINETIC THEORY 

Before discussing the details of LB methods, let us 
briefly review the theoretical background of kinetic the- 
ory in a relativistic context. The starting point is the 
Boltzmann equation for the single particle distribution 
/ = /(^jP ), with the relativistic analogue of the 
Bhatnagar-Gross-Krook collision term[4ll|: 



/ = - 



^ (/-/?), (i) 



Here a; M and p a are the position and momentum 4- 
vectors, respectively, tr is the single relaxation time, and 
/j q denotes the equilibrium distribution function. Here 
and in the following, we work in units where the speed 
of light c, the Boltzmann constant ks and Planck's con- 
stant h have been set to unity, c = fcs = ft = 1. In the 
ultrarelativistic case, this may be taken in the form of 
the Jiittner distribution function 
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where Z~ x parametrizes the number of degrees of free- 
dom, u M is the macroscopic (fluid) 4-velocity, T is the 
local temperature, V p denotes the (geometric) covariant 
derivative, and T^ v are the Christoffel symbols that are 
given by derivatives of the underlying metric tensor g^ v . 
The connection to fluid dynamics is realized by intro- 



ducing the energy- momentum (energy-stress) tensor T^", 

T^(t,x) = J d XP a p f(t,x,p) (3) 
d 4 p 



(27T) 



rS{p^ - m 2 )2H(p°)p a pPf(t,x,p) , 



where to is the particle mass and H the Heavyside step 
function. Note that we have introduced the notation 
f(t,x,p) = f(x^,p a ). The equations of motion obeyed 
by the tensor T^ v emerge, after a little algebra, upon in- 
tegrating Eq. ([TJ with respect to four-momentum degrees 
of freedom, J d\p v ', 



dxP^P" (/ /cq) = - — (T^ - T%) . 

TR TR 

The equilibrium energy-momentum tensor is read- 
ily computed using the Jiittner distribution, and reads as 
follows: 



cq 



(e + P)u»u u - Pg" 



where the energy density e and pressure P are functions 
of the temperature and the number of degrees of freedom 
Z (here arbitrarily set to one). The full energy momen- 
tum tensor may then be written as T^ v = T££ + W v , 
where the second term collects non-equilibrium contribu- 
tions. 

Still needed is the choice of the rest-frame of the heat 
bath with respect to which the fluid velocity is defined. 
In the ultrarelativistic limit, the canonical choice is the 
so-called Landau-Lifshitz condition, whereby u^T^ = 
eu v , so that Eq. reads simply as: 

expressing the (covariant) conservation of energy and mo- 
mentum. Sufficiently close to equilibrium, where gradi- 
ents are small, the form of can be calculated byin- 
tegrating Eq. (jl) with respect to / d X p^p a (c.f. [^El). 
In the ultrarelativistic case, when particle masses can be 
neglected (to = 0) one finds: 
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where 

A <« B I3> = 



and A MI/ = g^ v — u^u" . Performing a non-relativistic 
limit of V M T' 11/ one recovers the Navier-Stokes equations, 
with the following dynamic viscosity coefficient 



P 



V = TR- 



(i 



Therefore, Eq. ([I} reproduces the equations of fluid dy- 
namics in the continuum, on condition that tr = 6-T -1 , 
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where s denotes the entropy density and gradients must 
be small enough that a fluid dynamics description makes 
sense at all. A few remarks are in order: by con- 
struction, kinetic theory achieves a remarkable disen- 
tangling between non-linearity and non-locality, which 
proves beneficial for both theoretical and computational 
purposes. Indeed, in the hydrodynamic formulation, any 
generic quantity, including the flow velocity, is trans- 
ported along space-time changing trajectories, defined 
by the flow velocity itself, thereby giving rise to terms 
which are non-local and non-linear at a time. In a turbu- 
lent flow, such trajectories become typically fairly com- 
plicated, thus opening potential exposures to numerical 
inaccuracies and instabilities. In kinetic theory, on the 
contrary, information is always transported along con- 
stant characteristics, dx = vdt, since the velocity v (vec- 
tor notation relaxed for simplicity) does not depend on 
space-time coordinates. Thus, no matter how wild the 
space-time dependence of the fluid flow, the streaming 
operator is linear in the Boltzmann distribution function, 
and the information always travels along straight lines. 
The price to pay for this major advantage is the need of 
three extra-dimensions in velocity space. However, ve- 
locity space lends itself to very economic discretizations, 
typically 0(3 d ) discrete velocities in d spatial dimensions 
space, which make the tradeoff between non-linearity and 
over-dimensionality an excellent bargain, and is one of 
the key assets of the LB formulation. 

One may then wonder where the non-linearity has dis- 
appeared in the kinetic formulation. It turns out that it 
is concealed in the local equilibrium distribution, which 
is a non-linear function of the local hydrodynamic vari- 
ables, see expression @. Furthermore, note that since 
collisions are zero-ranged in the Boltzmann picture, the 
corresponding collision term is completely local, as antic- 
ipated. This lies at the root of the excellent amenability 
of LB to parallel computing, another major practical as- 
set of the technique altogether [32j . 

Having clarified the main philosophy of the kinetic 
pathway to fluid dynamics, one must come down to spe- 
cific details. The main question is: how sparse can the 
sampling in momentum space be made? 

The target criterion is that V^T^ = must emerge 
as a continuum limit of the lattice analogue of Eq. (TTJ). 
Like in the non-relativistic framework, this sets a spe- 
cific demand on the symmetry of the lattice tensors, as 
detailed in a sequel to this work. In particular, second 
order tensorial identities of the form 

u» J ' dxP^P^fit^x^p) = eu v . (5) 
have to be reproduced exactly in the lattice formulation. 

IV. DETOUR: THE NON-RELATIVISTIC 
LATTICE BOLTZMANN FORMULATION 

At this point, it is instructive to review the setup of 
the LB scheme in the non-relativistic context. There, the 



equilibrium distribution function for an ideal gas is the 
Maxwell distribution, 

where c s — y/T/m is the sound speed, the velocities v 
(alternatively denoted as v l with i running on spatial 
coordinates) take the role of the momenta p, and u is 
the macroscopic velocity. By introducing a scaled tem- 
perature 9 = ijr, To being a reference temperature, and 
rescaling the velocities with the speed of sound, the ve- 
locity moments (the conditions to recover fluid dynamics 
corresponding to Eq. ((SJ), take the form [22j 

J d 3 vf(t,x,v) = p, J d 3 vvf(t,x,v) = pu, (6) 

J d 3 vv 2 f(t, x, v) = 2pe int + pu 2 , 

where p is the mass density and £i n t the internal en- 
ergy density. Note that the local Maxwellian corresponds 
to the generating functional of the Hermite polynomials 

H n (v), 

e - ( ^ )2/2 = e-" 2/2 £i/„(«)^, (7) 

n 

in one-dimension. This is readily generalized to three 
dimensions, thanks to the factorizability of the pre- 
factor e~ v "/ 2 into the respective components. The actual 
Maxwellian distribution can then be approximated as 

^ TV. 

n=0 

which is valid in the sense of mean convergence, if 

J dve" 2 / 2 e-l v - u l 2 / 2e , 

exists, or equivalently 9 < 2 (c.f. [H, The 
same representation can be applied to the full distribu- 
tion f(t,x,v), which reads (after rescaling v — > v9 x / 2 , 
u^u9 1 ' 2 ): 

f{t,x,v) = e-v 2 ' 2 hm f;<^M^-«-(v), (8) 

where W£ ■ l " are tensor Hermite polynomials (c.f. 23]). 
In practice, / will be approximated by truncating the 
above sum at finite (small) N. 

Integrals of the form J e~ v / 2 P(v)dv, where P is a 
polynomial of degree 2N or less, can then be calcu- 
lated exactly as a sum over the roots of Hn(v) ("Gauss- 
Hermite quadrature"). Therefore, the roots v = v m , 
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D2Q9 D3Q19 



FIG. 1. Typical lattice configurations D2Q9 (9 velocities in 
2 dimensions) and D3Q19 (19 velocities in 3 dimensions) for 
the lattice Boltzmann model. The arrows denote the discrete 
unit vectors set. 

m = 1 , . . . N represent the ideal choice for the discretiza- 
tion/sampling of velocity space. To guarantee that the 
Boltzmann equation reproduces the non-relativistic fluid 
dynamics equations, the velocity moments must be rep- 
resented exactly, which implies the necessary condition 
N > 2. For many applications, specifically those not 
dealing with strong thermal and compressible phenom- 
ena, this is also sufficient to a second order numerical ac- 
curacy. The end result is a discrete (lattice) Boltzmann 
equation of the form 

A(f + At J x + c i At)-/ i (i,x) = -fi(/i-./H, (9) 

where fl = ^ with tr the single relaxation time. The 
discrete velocities c, run over a lattice with sufficient sym- 
metry to guarantee mass, momentum and momentum- 
flux conservation ([6]) , so as to recover the exact form 
of the Navier-Stokes equations. Typical lattices fulfill- 
ing the above constraints are the D2Q9 (nine velocities 
in two dimensions) and D3Q19 (nineteen speeds in three 
dimensions), see Fig. [TJ 

Note that, owing to the Cartesian formulation, these 
lattice configurations are space-filling. This is crucial 
to ensure the so called light-cone condition dxj = Cidt, 
i.e. the discrete populations hop from site to site in 
fully synchronous mode (grid-bound dynamics), a fea- 
ture which turns out to be critical to the computational 
efficiency of the LB scheme. Upon performing a stan- 
dard Chapman-Enskog asympt otic expansion, the lat- 
tice Boltzmann scheme ( |28l - [30| ) is shown to recover the 
hydrodynamic equations of a quasi-incompressible fluid 
with kinematic viscosity (in lattice units At = Ax = 1) 
v = c 2 s (tr — 1/2), c s being the lattice sound speed, 
typically l/y/3. To be noted, the factor —1/2 at the 
right hand side, which stems from a second order Tay- 
lor expansion of the discrete streaming LBE operator. 
This term, known as 'propagation viscosity', contributes 
a negative viscosity and permits to achieve very small 
viscosities v -C 1 with unit time-steps, by simply choos- 
ing tr = 1/2 + v/(? s . This property is crucial to access 
low- viscous, turbulent regimes while still preserving an 
efficient time-marching procedure. 



V. A FULLY RELATIVISTIC LB ALGORITHM 

The lattice formulation of relativistic kinetic theory 
poses a few genuinely new challenges, primarily the fact 
that the energy E is no longer a simple quadratic function 
of the velocity (momentum) E — y/ m 2 + p 2 . This ba- 
sic feature reflects in the non-separability of the Jiittner 
distribution along the three components of the momen- 
tum p, and forbids a simple transcription of the Hermite 
procedure described above for the case of non-relativistic 
fluids. 

This is the reason why the only existing relativistic 
LBE version available to date is based on a top-down pro- 
cedure, i.e. design lattice equilibria with free Lagrangian 
parameters, which are then adjusted in such a way as 
to match the five basic conservation laws, number den- 
sity and energy-momentum 4-vector. Full details can be 
found in the original references [HI, Ho[. The scheme 
was validated for two different relativistic applications, 
Id quark-gluon plasmas, 3d supernova explosions, and 
graphene [42j, showing excellent performance on all of 
them. However, inherent to the moment-matching pro- 
cedure, is the question of realizability, i.e. the existence 
of an underlying microscopic model or at least an equiv- 
alent analogue in continuum kinetic theory. Moreover, 
the top-down procedure only works in the case of Carte- 
sian coordinates. In view of general relativistic applica- 
tions, involving generic coordinate systems, it is highly 
desirable to explore the viability of the relativistic LB 
procedure beyond the Cartesian realm (incidentally, this 
would prove useful also for non-relativistic applications 
in general coordinates). 

For the ultra-relativistic case, where particle masses 
can be neglected, the equilibrium distribution reads as 

e -p-u/T = e -|p|u°/T+p-u/T ^ £ 1Q \ 

which does not allow an ansatz such as Eq. ([7} because 
the pre-factor corresponds to an unfactorizable square 
root dependence. This alone prevents a straightforward 
use of Cartesian coordinates. 

Rather, the equilibrium distribution function suggests 
the following expansion: 

P . u/T = M u°/(T e) sp(P_Y( \p\f\ * (")" 

^\\p\J \T e) («°)»n!' 

(11) 

which involves unit vectors v = t^t rather than momenta 



p and powers of jFg-, that go together with the expo- 
nent and where again the scaled temperature 9 — ^- was 
introduced. 

This leads to the following ansatz for the relativistic 
distribution function /: 

/(*,*,„) =e-^^J2 P S.i n (y)^ An (t,x,p°/T ) , 

n 

(12) 



5 



where p° = |p| and vector polynomials P^ n \ (v) which 
are orthogonal with respect to the angular integral J . 
Their properties are listed in appendix [AJ 

The ansatz (fT2|) should constitute a valid approxima- 
tion to f(t,x,p) in the sense of mean convergence, i.e. 
provided that the integral 



exists, which together with Eq. (fT0| implies 9 < 
2 (u° — |u|). For highly relativistic fluid flows such as 
those with Lorentz factor of 7 ~ fO, 9 = 0.1 would ensure 
convergence. However, in practical applications, such as 
the dynamics of heavy-ion collisions, such high fluid flow 
velocities only occur for the small temperature region. 
Hence, setting To to the maximum temperature encoun- 
tered in the problem was found to give acceptable results. 

Denoting the scaled momentum as p = p° /To, it proves 
convenient to further expand the coefficients oW in terms 
of generalized Laguerre polynomials , so that the 
complete ansatz for / is given by 



f(t,x,p)=e-P 



N p -1 Nv-1 

E E 

k=0 n=0 



where the choices a = 2, 3 will be most relevant, so that 
a is restricted to integer numbers hereafter. For conve- 
nience, the main properties of Laguerre polynomials are 
reported in appendix iBl Using orthogonality, the coeffi- 
cients cS nk ^ up to second order read as follows: 



(nk) 



(m + a) ! 

3m! 
[m + a)\ 
15m! 
2(m -I- a)! 



dpp a 



dpp a 



4ir J m 


(P) 


= a ( 0m 


[t,x), 


— fP^L^ 


(P) 


(lm 


(t,x) , 


^ fp (2) L (a 
4tt J « m 


(P) 


(2m 
= % 


M14) 



where we have used the fact that a? may be taken to 
be traceless and symmetric. 

The ideal choice for the discretization of velocity space 
is once again given by the requirement that the velocity 
moments ([5]) be represented exactly. Use of non Carte- 
sian coordinates, however, implies that the associated lat- 
tice structure is no longer space-filling in general. 

For a polynomial P(p) of degree less than 2N P , the 
integral over p can be recast into an exact sum 



[ dppe-fp(p)= jr u$P(p 
J »=o 



if the nodes po, ■ ■ - Pn p -i are given by L^(pj) — and 
the weights are given by 



p _ (JVp + a)! 



Pi 



Npl (N p + 1)2 



-(«) 

J N P + 



l(Pi) 



Note that the requirement ([5]) implies a polynomial P{p) 
of degree 3 — a + N p , implying that N p > 3 — a is a 
necessary condition to represent this polynomial exactly. 
The choice of a may depend on the problem, but a = 3 
is particularly convenient here because it minimizes N p . 
Unless one is interested in considering finite chemical po- 
tential, it is therefore convenient to choose a = 3, i.e. 
Np = 1, as we shall do from now on. 

Coming to the angular dependence, for polynomials P 
of sm<j) and cos</> of degree less than 2N$, integration 
over the polar angle </> can be expressed as 

/ d<f>P(cf>) = Y, <»tP(<M > 

J -"* 1=0 

where <f>i = an d ujf = A polynomial p( n \v l ) 
will only involve powers of cos cj>, sin <f> up to n, so an 
exact representation of Eq. ([5]) requires N v > 3, hence 
> 3. 

Finally, for the integration w.r.t to the cosine of the 
polar angle cos(# p ) = £, there is a special symmetry that 
one can exploit. Namely, should the integrand contain a 
p x /p° = cos^yl — £ 2 , then the ^-integration will only 
give a non-zero value if another cos <f> is present in the 
integrand. The only way this can happen is through 
another factor p x /p°, which means that cos 2 <fi (l — £ 2 ) 
must be present in the integrand. A similar argument 
may be given for p y . Hence, the integral over the po- 
lar angle is always of the form Jj^ 1 d£P(£), with P(£) a 
polynomial of degree 2 + N v or less. This is represented 
accurately as 



iv t -i 

J=0 



(15) 



for polynomial degrees less than 2N^, £j being the roots 
of the Legendre Polynomial P/v £ (x) and 



■u'i 



(l-3)(^ t &) 



(16) 



the corresponding weights. 

The requirement ([5]) again implies N% > 3. 

To summarize, the requirements (O are fulfilled for the 
ansatz (1131) if one uses a discretized set of momenta 



Piji 



T pt 



-(i) 



/ f 

cos 4>i , /l 



sm^.,/1 
6 



/ 



the roots of Pw, and = 



with pi the roots of L^* +1 , Xj 
j^j- , with JV P , iVg , greater or equal to 3 and 3 < iV„ < 
min(A^,,iVj) (see Fig.[5]for illustration). 
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as follows: 



FIG. 2. Example of a lattice configuration for the present rel- 
ativistic model. The arrows denote the discrete spatial com- 
ponents of p'\ 



This implies a minimum number of 27 discrete 
"speeds" j/ 1 ^;, quite comparable with the number of dis- 
crete speeds commonly used in non-relativistic LB theo- 
ries. 

With the momentum space thus discretized, we can 
move ahead and set up a concrete, fully relativistic lattice 
Boltzmann algorithm. 

Before doing so, however, it is instructive to study how 
to extract the fluid velocity and local temperature distri- 
bution from a given distribution function / = fiji(t,x) 
discretized via (11; 



A. Energy-Momentum Conservation 



Using Eq. (|13p , the energy momentum tensor becomes 




,(20) 



(17) 



The requirement ([5]) implies that should be a future- 
pointing eigenvector of the energy momentum tensor, or 
u^T^ — eu" . Thus, one has to calculate the eigensys- 
tem of and identify w M as the (only) future-pointing 
eigenvector, with eigenvalue e. Using existing numerical 
packages, this can be done in a rather efficient way. Once 
e is known, the temperature is calculated from e = 
Since we neglected particle masses, the equation of state 
is always that of an ideal gas, or P = |. Non-ideal equa- 
tions of state will be considered in a follow-up work. 



B. The equilibrium distribution function 

While the full equilibrium distribution function is given 
by Eq.©, the LB algorithm approximates all /'s by the 
ansatz (fT5|) . Hence, for consistency, / oq is also expanded 



Hit, x, p) = e-P (t, x) + Pia^l (i, x) 



+ -Pij a ijjq(t> x ) 



(18) 



where we neglect higher order terms. The coefficients are 
readily evaluated to be 



a 



(oo) = 

cq 
(10) 



1 + -u' 



a 



l.cq 
(20) 



(30) 



Id ( u v_u 2 ^ 



log (1 + 2u 2 + 2|u|u°) )e 
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With these expressions at hand, we are now ready to 
construct an operational algorithm. 



VI. THE FULLY RELATIVISTIC LB 
ALGORITHM 

First of all, a lattice version of the Boltzmann equa- 
tion ([T]) needs to be established. For this purpose, we 
define three steps: i) Streaming in configuration space 
(cc-move), ii) Streaming in momentum space (p-move), 
iii) Collisional relaxation. These read as follows (discrete 
indices are relaxed for notational simplicity): 

f'(t,x i ,p a ) = f(t + 8t,x i + ^5t, J A , (19) 

f"{t,x\ P a ) = f( tl x\p a ) + 5tr^^d^f(t,x\ P a ) , 



f'(t,x\p a ) = f(t,x\p a ) 



p^u^St 

P°TR 



if" - / oq ) 



Here, we use the notation f(t, x l ,p a ) = f(t, x,p) to make 
the vector components explicit. 

The first step consists in the streaming of the distribu- 
tion functions according to the discrete momenta. The 
second one, is the implementation of the external forces 
due to the curvature of the space-time, and the third 
one is the collision step, expressed in terms of relaxation 
towards the local equilibrium. 

The relativistic lattice Boltzmann algorithm is given 
by the following sequence of steps: 

1. Initialization: 

At the initial step, /' needs to be known at grid 
sites x l and discretized momenta p a . For this, we 
introduce the initial conditions of the specific prob- 
lem and the initial distributions is typically speci- 
fied as the local equilibria associated with the initial 
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hydrodynamic fields f(t = to,x l ,p a ) = f eq (t,x,p), 
where the equilibrium distribution is given by the 
expression (fl~8|) . 

2. x-move: 

Calculate the new / from 

f(t,x\p a ) = f (t~6t,x*-v l 5t,p a ). 

In non-relativistic lattice Boltzmann methods, this 
step is known as streaming, each populations moves 
to the site pointed by its corresponding discrete 
speed. 

3. p-move: 

Compute the change in the distribution function 
because of external/internal forces. In princi- 
ple, a force term implies a change in velocity, 
thereby jeopardizing the discrete structure of ve- 
locity space. However, this can be preserved by 
moving the information according to the streaming 
step given above (with constant speeds/momenta) 
and then correcting the distribution form with an 
appropriate source term. The latter is identified 
by representing the derivative in momentum space 
also as a polynomial expansion: 

N v -1 
n=0 

(20) 

with F^ = T^^f. 

The unknown source coefficients s'"' can be com- 
puted by inverting (|2"0"j) and using integration by 
parts. Note that in general these coefficients can 
be expressed as sums over the coefficients a("°) in 
Eq. PH. 

With the source term available, we evaluate /" as 
f"{t,x\p a ) = f{t,x\p a ) 

n=0 

This step accounts for the geometric forces in the 
given space-time. 

4. Equilibration: 

In order to perform collisional relaxation, local dis- 
crete equilibria must be constructed first. To this 
purpose, we calculate the energy momentum tensor 
T^ v corresponding to /" from Eqs. (HI) and (fH|) . 
We then compute the values of the fluid 4-velocity 
it M and energy density e by calculating the future- 
pointing eigenvector of T^ v . The local temperature 
is obtained by the equation of state. Using T, it^, 
we calculate / oq from Eq. (fT5)) . and the collision 
coefficient £1 = p fl u ti /(p°Tn). 



5. Collision: 

Calculate the post-collisional state /' from the 
known values of /", / cq , and Jl, according to 

f(t,x\ P a ) = f'(t, x\ P a )(i~nst) 

+ n8tf c «(t,x\p a ). (21) 

6. Cycle through 2 — 5 for each time-step until com- 
pletion of the time evolution. 

A remark concerning the x-move step is in order. Since 
the unit vectors v l — p l /\p\ discretize the unit sphere (see 
Fig. [5]), the displaced positions x % — v l St will typically not 
correspond to a neighbouring grid site, unless very par- 
ticular geometries (e.g. a hexagonal lattice) are chosen. 
This breaks the light-cone rule discussed previously. 

At least two ways out of this problem can be envisaged. 

The first is to acknowledge the fact that spatial and 
momentum discretization can no longer be kept in syn- 
chrony, going back to the original Boltzmann equation ([1]) 
and discretize space derivatives in a flux-conserving way, 
according to one's favored finite- volume/difference prac- 
tice. This is similar to the non-relativistic LB method 
on so-called unstructured meshes, wherein powerful fea- 
tures of modern finite-volume techniques are imported 
within the LB framework. For more details about this 
technique, we refer the reader to the literature, e.g. 
Refs. [43|, |44[ . The second method, the one we will adopt 
in the following, consists of transferring off-grid distribu- 
tions into grid locations through (bilinear) interpolation. 

Either way, it is clear that none of the two methods 
above can match the simplicity, hence computational ef- 
ficiency, of the space-filling Cartesian formulation. In 
particular, they cannot preserve the exact nature of the 
streaming step in light-cone form. This limitation ap- 
pears to be inherent to the non-separability of the rel- 
ativistic Jiittner equilibria along the three spatial coor- 
dinates. In this respect, a complete transfer of the key 
assets of the non-relativistic scheme, does not appear to 
be feasible. This is the price to pay for a fully rela- 
tivistic lattice formulation. However, since information 
is still traveling along constant streamlines in configura- 
tion space, this does not prevent the fully relativistic LB 
algorithm from delivering competitive performance, as 
we shall demonstrate in the next section, where concrete 
test-case simulation are presented. 



VII. APPLICATION: THE BOOST-INVARIANT 
QUARK-GLUON PLASMA 

As a first application, let us consider the Milne space- 
time generated by the coordinate transformation r = 
y/t 2 — z 2 , Y = tanh _1 (z/i). In these coordinates, the 
metric takes the form 

= di&g(g TT ,g xx ,g yy ,g Y Y) = diag(l, -1, -1, -r 2 ) , 
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where hereafter the (+,—,—,—) sign convention is as- 
sumed. In this metric, the non-vanishing Christoffel sym- 
bols are: 



r Y -T Y - - 
L tT - l Yt - T 



Tyy=T. 



This implies a non-vanishing covariant fluid divergence 
even for a fluid a rest, u M = (1,0), i.e. 



r 



i 



The reason for this is that the Milne space-time is ex- 
panding in one dimension, so that a system at rest ex- 
periences the 'stretching' of space-time. This is a nice 
feature, because it naturally mimics the expansion of 
the system, following a heavy-ion collision in Minkowski 
space (see Ref. [24[ , Section 5B for details) . In Minkowski 
space-time, a general solution to p^d^f — is e.g. 
/ = f(p±,tp — xp f ). In Milne-coordinates that would 
correspond to 

/ = f(p±, r 2 p Y , cosh Y (rp± — xj_p T ) — rp x± sinh Y) . 

(22) 

This can be further simplified by considering only the 
evolution at mid-rapidity, Y ~ 0. It is readily checked 
that, under such condition, the action of the derivative 

dr is exactly cancelled by the dy derivative of the last 
term in (|22p. One may thus neglect both, so that the 
Boltzmann equation at mid-rapidity simplifies to 



p T d T +p±_ ■ d± 



-c[f], r~o, 

(23) 

where p T is treated as an independent variable (e.g. 

dy p r — 0). The corresponding discrete version of the 
lattice Boltzmann equation takes the following form: 



/ ( T + 5T,X l + Z-6T, P a 

P T 



f'{r,x\p a ) 



(24) 



f"(r,x\p a ) = f(r,x\p a ) + —p Y d^f, 

T 

f , (r,x\p a )=f"(r,x i ,p a ) 



P t tr 



if" - D 



A. Warmup: O+ldimensions in Milne space-time 

The simplest practical example is given by considering 
a system that is homogeneous in the Milne coordinates, 
e -g- / — f( T >P a ) = f( T iP,£i)- Because there is no space- 
dependence left, one may use a simplified version of the 
discretization ansatz (|13l) . 



/(r,p«) = e- 4 ^0k(r)P k (0 , 



(25) 



where p z — p t, such that p T — yj p\ -I- (p y r) becomes 

p T = |p| = \Jp\ +p1 and £ — arccos^. With / dis- 
cretized this way, one readily inverts to obtain the coef- 
ficients a m i. For instance 



2/ + 1 
12 



Nf-1 



3=0 



(26) 



where the second line is the discretized version and we 
recall that £j are the roots of Pn 6 and the weights u>* 
were specified above. Starting with an initial equilibrium 
distribution function 

/(T ,p,&)c< e -^, 

and discretizing /(ro) on a £ grid with N% points (p = 4), 
the algorithm reads as follows. 
First set /' = /. Then, set 

where the coefficients are calculated to be 

s° = ~j dp? J da{i+e)f, 



S 2 = - 



24 



Calculate /" from 

f"(-n+5T,p,ti) = /(T ,p,£)+^Ve-P (S° + S 2 P 2 (0) 

T + Or 

Next, calculate the equilibrium temperature via the 
energy-momentum tensor corresponding to /" and ob- 
tain f (r + St) from Eq.([2l]). 

The above steps are cycled in time over the prescribed 
time-span of the simulation. 



B. Comparison with exact results 

Let us now compare the results of the above lattice 
Boltzmann algorithm in Milne-spacetime against known 
exact results. First, let us consider the case in which the 
relaxation time is so large that the collision term plays 
a negligible role. The solution to the Boltzmann equa- 
tion must therefore be very close to the free streaming 
solution, 



(r,p,e)=e^V^/Q. 



/free— stream 
2 

where r = ^ — 1 and Q is the initial temperature scale. 
From this, the temperature can be obtained from the 
energy-momentum tensor as 

1/4 

I / T-f 1 nrntnn , / t z /t£ — I \ 

T(t) = Q 



k=0 



Fq arctan \/t 2 /tq — 1 



\A7r 2 -i 
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As opposed to the free-streaming case, let us now con- 
sider the opposite extreme of small viscosity, i.e. very 
fast relaxation to the local equilibrium. In this case, the 
'exact' temperature evolution is given by fluid dynam- 
ics. More specifically, denoting Ty — P = <£, the energy 
density and $ fulfill the coupled equations [24( 



d T e = 
<3 r $ = 



P 



T 



TV 



Arj 
3t t t 



4$ 
37 



Ai 



(27) 



where t„ is the relaxation time and Ai is a self-coupling 
parameter that can be calculated [25j , presumably fairly 
easily for the BGK collision kernel used here. 

Since both are second-order corrections to hydrody- 
namics, their determination is left for future work. Here, 
we simply set Ai = / 6 '' m and vary tv between tv 

2(2^*2)2 ^ 



i s cmiu — &^ (the weak and strong coupling 
limit, respectively, [21|). The above set of hydrodynamic 
equations display simple analytic solutions for the case of 
vanishing viscosity (ideal fluid) and first order gradient 
expansion (Navier-Stokes) . These are [13, [H[ given by 



T(t).T (^ 1/3 



i + ^(V(^ 2/3 

3st T \ V T 



(28) 



The full set of equations (|27p is second-order in gra- 
dients and thus causal for tv larger than some critical 
value. For general values of T n ,\i, it has to be solved 
numerically. 

In FigJ21 we show a comparison between the lattice 
Boltzmann algorithm, against various 'exact' results for 
the case of rj/s = 0.5 and T (t = lfm/c) = 0.5 GeV. 

The second-order set of hydrodynamic equations was 
solved using forward time differencing St8 t X(t) = 
(X(t + St) - X(t)). A time step of St = 0.01r was 
required to reach a stable continuum result. Conversely, 
for the lattice Boltzmann algorithm, typically the result 
is stable for St < 0.2to, nearly 20 times larger than 
the fluid dynamics requirement and only about 5 times 
smaller than standard LB schems in cartesian geome- 
tries. Based on the general arguments discussed earlier 
on in this paper, we interpret this encouraging outcome 
as the beneficial effect of moving information along con- 
stant streamlines. 

As can be seen from this figure, the lattice Boltzmann 
algorithm does track the 2 nd order viscous fluid dynam- 
ics result from early times onwards. (The Navier-Stokes 
result has a different initial condition and hence the other 
results are not expected to track this curve). The num- 
bers in Fig. [3] were chosen such that the initial tempera- 
ture corresponds to the maximum temperature expected 
for heavy- ion collisions at -y/s = 5.5 TeV at the Large 
Hadron Collider. Moreover, T — 0.15 GeV is the temper- 
ature where a freeze-out to hadrons is expected. Hence, 
for this temperature region, the lattice Boltzmann algo- 
rithm with N$ = 5 accurately reproduces the 2 nd order 
viscous fluid dynamics results. 



Temperature 



0.5 



0.125 



T=0.15GeV 


■ T)/s = oo 

r|/s = 0.5 (l sl order) 
• q/s = 0.5 (LB) 
q/s = 

■ ■ ■ r\/s = 0.5 (2 nd order) 
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X [fm/c] 



100 



Scaled Temperature 



1.3 



1.2 



■ T|/S 




r|/s 


= 0.5 (1 st order) 


■ ■ ■ T|/S 


= 0.5 (2 nd order) 


• T|/S 


= 0.5 (LB) 



T=«).15 GeV. 




100 



FIG. 3. Temperature evolution of a system starting at tq — 1 
fm/c and temperature To = 0.5 GeV. Results shown are 1 st 
order viscous fluid dynamics (Navier-Stokes) (initially not in 
equilibrium), 2 nd order viscous fluid dynamics (two values of 
Ttt, see text), and the lattice Boltzmann result (LB) for = 
5, all for r)/s — 0.5. As a reference, the free streaming result 
(r]/s — > oo) and the result for ideal fluid dynamics (rj/s — 0) 
are also shown (top panel). Bottom panel: results are divided 
by the ideal fluid dynamic result to highlight differences. The 
LB result is seen to converge to the 2 nd order viscous fluid 
dynamics. 



The dependence of the lattice Boltzmann result on the 
chosen discretization is shown in FigHJ There, one can 
see that all cases (even JVj = 3) reproduce the viscous 
fluid result, and for finer discretization > 4, the LB re- 
sults arc indistinguishable to the naked eye. The viscosity 
dependence is also shown in FigJU The lattice Boltzmann 
algorithm tracks the fluid dynamic result for viscosities 
up to rj/s ~ 1.0. For smaller viscosities, rj/s < 0.5, LB 
undershoots the fluid dynamics result. However, decreas- 
ing St by a factor two, brings the LB back in line with 
the fluid dynamic result. 
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N , dependence 

T\/s = 0.5 




T|/s dependence 




100 



Temperature profiles 

b=7 fm collision, r|/s=0.08 
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FIG. 4. Top: dependence of temperature evolution on dis- 
cretization (for r\/s — 0.5). Bottom: dependence of result 
on value of shear viscosity coefficient compared to 2 nd order 
viscous fluid dynamics (for TVj = 5). Viscous fluid results 
from (|27|l for two values of tv, see text. For the lowest vis- 
cosity shown (r]/s = 0.2) we highlight the effect of numerical 
viscosity by choosing two different values of St. 



C. 2+1 dimension in Milne space-time 

Next we consider the case in which the transverse dy- 
namics is also taken account. Inclusion of the transverse 
space dynamics requires the solution of Eq. (|24[) . using 
the full discretization (fT5|) . For convenience, we choose a 
square lattice for the grid in x±. Choosing furthermore 
St = 8x, we use bilinear interpolation to obtain /' at 
points that lie in-between lattice sites. The change in 
momenta p is calculated similar to Eq. (p?0"|) . 

The LB solver is applied to simulate the evolution of 
the medium created in Au + Au collisions at top RHIC 
energies (y/s — 200 GeV). For smooth initial conditions, 
the results may then be compared to the fluid dynamics 
solution, given by the code VH2+1 This code has 

been cross-tested against several other codes and is gen- 
erally credited for producing reliable results for smooth 
initial conditions. 

The initial conditions are generated at initial time r = 



FIG. 5. Top: temperature evolution in viscous hydrody- 
namics ('hydro') versus lattice Boltzmann equation ('LB') for 
rj/s = 0.08. As can be seen, the temperature evolution in 
the lattice Boltzmann approach for Sx = 0.2 fm is reasonably 
close to the 'exact' hydrodynamic result. Bottom: evolution 
of velocity u x /u° in viscous hydrodynamics versus LB. Even 
high velocities up to 80 percent of the speed of light are well 
represented, but at later times an instability develops in the 
low temperature region (a grey triangle marks the position 
of T = 0.15 GeV). Note that the discrepancy at larger x lies 
exclusively in the region T < 0.15 GeV (compare top plot). 



1 fm/c from a Glauber model, with number of collisions 
scaling (c.f.0) on a 69 x 69, 139 x 139 or 279 x 279 grid 
with a lattice spacing of Sx = 0.4 fm, Sx — 0.2 fm or 
Sx = 0.1 fm, respectively. The maximum temperature 
at the center of the grid is T max = 0.37 GeV for central 
collisions (impact parameter 6 = fm). 

In Fig. lVIICl the temperature evolution from LB with 
= 4, — 4 for Sx = 0.2 fm is compared to the 
VH2+1 solution. As can be seen from this figure, the LB 
code reproduces the VH2+1 result rather accurately for 
temperatures T > 0.15 GeV, the relevant temperature 
regime for the fluid medium. At later times r > 5 fm/c, 
when all fluid cells have cooled below a temperature of 
0.15 GeV, numerical instabilities develop at the outer 
edges, for the discretization used N% = 4, = 4. We 
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have checked that the remaining discrepancy between 
fluid dynamics and LB in the temperature evolution close 
to x ~ can be cured by using a smaller lattice spacing 
Sx = 0.1 fin. 



Also shown in Fig. I VII CI is the evolution of the veloc- 
ity One can see that LB tracks the VH2+1 result 
closely for the high temperature region (there are clear 
deviations at low temperatures T < 0.15 GeV, compare 
left plot in Fig. I VII Cj) . Even high velocities seem to 
be well represented, but the algorithm does not handle 
correctly velocities u x /uq > 0.8. Presumably, this dis- 
crepancy can be cured by including higher order terms 
in Eq. fTE)) . However, it should be pointed out that for 
the high temperature region T > 0.15, the velocity distri- 
butions from fluid dynamics are accurately reproduced. 



VIII. CONCLUSIONS 



Summarizing, we have developed a new scheme based 
on the lattice-Boltzmann method to model relativistic 
fluid dynamics in general spacetime. The main advan- 
tage of our scheme, as compared with previous relativistic 
lattice Boltzmann models [1^, H(| , rests mainly with its 
ability describe the dynamics of ultra-relativistic systems 
in general space-time geometries. The present model dif- 
fers from typical lattice Boltzmann schemes mostly in the 
streaming step, which, because of the spherical shape of 
the discrete momenta, is no longer space-filling. Instead, 
multi-linear interpolation is used to represent the distri- 
bution functions in the second-nearest neighbours of each 
cell on the lattice. This interpolation breaks the exact na- 
ture of the standard LB streaming operator. However, at 
variance with hydrodynamic formulations, it still moves 
information along constant streamlines, thereby permit- 
ting to march in larger time-steps than hydrodynamic 
codes. Our scheme has been validated through simula- 
tions in quark-gluon plasma, yielding very satisfactory 
agreement with other computational methods based on 
a macroscopic description, at a lower computational cost 
(nearly two orders of magnitude faster than the VH2+1 
viscous hydro code [6| and still a factor 3 — 5 as compared 
to optimized ones |45|). 

Because of these favorable properties, we expect this 
new LB method to offer a new competitive entry for the 
computational study of large-scale complex relativistic 
fluids. 
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Appendix A: The p( n ' Legendre polynomials 

The vector polynomials -P/™^ (v) are constructed by 
requiring orthogonality with respect to the angular in- 
tegral J 4p . One finds Specifically, the polynomials in- 



volved are £ 


;iven by 


p(0) 


= 1, 


p (i) 

i 


= Vi , 


p(2) 
ij 


1 

= ViVj - -■ 


p(3) 


= ViVjVk - 



(Al) 



The orthogonality relations for the first few polynomials 
are found to be 

47T 

^pU)p(i) = 5 j± 
Air 1 j 3 ' 

(A2) 



Appendix B: Generalized Laguerre Polynomials 



The first few generalized Laguerre Polynomials are 



;iven by 




r(«) 








L 2 




r(«) 





(Bl) 



x 3 (a + 3)x 2 (a + 3)(a + 2)x 
~~6 + 2 2 
(a + 3)(a + 2)(q + l) 
+ 6 
The orthogonality relation is given by 

1 dx e~' x L^\x)L^\x) = T{n + C ; + 1) 6 nm 



(B2) 
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